home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpklib.zip / SSICO.FOR < prev    next >
Text File  |  1984-01-07  |  8KB  |  246 lines
  1.       SUBROUTINE SSICO(A,LDA,N,KPVT,RCOND,Z)
  2.       INTEGER LDA,N,KPVT(1)
  3.       REAL A(LDA,1),Z(1)
  4.       REAL RCOND
  5. C
  6. C     SSICO FACTORS A REAL SYMMETRIC MATRIX BY ELIMINATION WITH
  7. C     SYMMETRIC PIVOTING AND ESTIMATES THE CONDITION OF THE MATRIX.
  8. C
  9. C     IF  RCOND  IS NOT NEEDED, SSIFA IS SLIGHTLY FASTER.
  10. C     TO SOLVE  A*X = B , FOLLOW SSICO BY SSISL.
  11. C     TO COMPUTE  INVERSE(A)*C , FOLLOW SSICO BY SSISL.
  12. C     TO COMPUTE  INVERSE(A) , FOLLOW SSICO BY SSIDI.
  13. C     TO COMPUTE  DETERMINANT(A) , FOLLOW SSICO BY SSIDI.
  14. C     TO COMPUTE  INERTIA(A), FOLLOW SSICO BY SSIDI.
  15. C
  16. C     ON ENTRY
  17. C
  18. C        A       REAL(LDA, N)
  19. C                THE SYMMETRIC MATRIX TO BE FACTORED.
  20. C                ONLY THE DIAGONAL AND UPPER TRIANGLE ARE USED.
  21. C
  22. C        LDA     INTEGER
  23. C                THE LEADING DIMENSION OF THE ARRAY  A .
  24. C
  25. C        N       INTEGER
  26. C                THE ORDER OF THE MATRIX  A .
  27. C
  28. C     OUTPUT
  29. C
  30. C        A       A BLOCK DIAGONAL MATRIX AND THE MULTIPLIERS WHICH
  31. C                WERE USED TO OBTAIN IT.
  32. C                THE FACTORIZATION CAN BE WRITTEN  A = U*D*TRANS(U)
  33. C                WHERE  U  IS A PRODUCT OF PERMUTATION AND UNIT
  34. C                UPPER TRIANGULAR MATRICES , TRANS(U) IS THE
  35. C                TRANSPOSE OF  U , AND  D  IS BLOCK DIAGONAL
  36. C                WITH 1 BY 1 AND 2 BY 2 BLOCKS.
  37. C
  38. C        KPVT    INTEGER(N)
  39. C                AN INTEGER VECTOR OF PIVOT INDICES.
  40. C
  41. C        RCOND   REAL
  42. C                AN ESTIMATE OF THE RECIPROCAL CONDITION OF  A .
  43. C                FOR THE SYSTEM  A*X = B , RELATIVE PERTURBATIONS
  44. C                IN  A  AND  B  OF SIZE  EPSILON  MAY CAUSE
  45. C                RELATIVE PERTURBATIONS IN  X  OF SIZE  EPSILON/RCOND .
  46. C                IF  RCOND  IS SO SMALL THAT THE LOGICAL EXPRESSION
  47. C                           1.0 + RCOND .EQ. 1.0
  48. C                IS TRUE, THEN  A  MAY BE SINGULAR TO WORKING
  49. C                PRECISION.  IN PARTICULAR,  RCOND  IS ZERO  IF
  50. C                EXACT SINGULARITY IS DETECTED OR THE ESTIMATE
  51. C                UNDERFLOWS.
  52. C
  53. C        Z       REAL(N)
  54. C                A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT.
  55. C                IF  A  IS CLOSE TO A SINGULAR MATRIX, THEN  Z  IS
  56. C                AN APPROXIMATE NULL VECTOR IN THE SENSE THAT
  57. C                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
  58. C
  59. C     LINPACK. THIS VERSION DATED 08/14/78 .
  60. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  61. C
  62. C     SUBROUTINES AND FUNCTIONS
  63. C
  64. C     LINPACK SSIFA
  65. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  66. C     FORTRAN ABS,AMAX1,IABS,SIGN
  67. C
  68. C     INTERNAL VARIABLES
  69. C
  70.       REAL AK,AKM1,BK,BKM1,SDOT,DENOM,EK,T
  71.       REAL ANORM,S,SASUM,YNORM
  72.       INTEGER I,INFO,J,JM1,K,KP,KPS,KS
  73. C
  74. C
  75. C     FIND NORM OF A USING ONLY UPPER HALF
  76. C
  77.       DO 30 J = 1, N
  78.          Z(J) = SASUM(J,A(1,J),1)
  79.          JM1 = J - 1
  80.          IF (JM1 .LT. 1) GO TO 20
  81.          DO 10 I = 1, JM1
  82.             Z(I) = Z(I) + ABS(A(I,J))
  83.    10    CONTINUE
  84.    20    CONTINUE
  85.    30 CONTINUE
  86.       ANORM = 0.0E0
  87.       DO 40 J = 1, N
  88.          ANORM = AMAX1(ANORM,Z(J))
  89.    40 CONTINUE
  90. C
  91. C     FACTOR
  92. C
  93.       CALL SSIFA(A,LDA,N,KPVT,INFO)
  94. C
  95. C     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
  96. C     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
  97. C     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
  98. C     GROWTH IN THE ELEMENTS OF W  WHERE  U*D*W = E .
  99. C     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
  100. C
  101. C     SOLVE U*D*W = E
  102. C
  103.       EK = 1.0E0
  104.       DO 50 J = 1, N
  105.          Z(J) = 0.0E0
  106.    50 CONTINUE
  107.       K = N
  108.    60 IF (K .EQ. 0) GO TO 120
  109.          KS = 1
  110.          IF (KPVT(K) .LT. 0) KS = 2
  111.          KP = IABS(KPVT(K))
  112.          KPS = K + 1 - KS
  113.          IF (KP .EQ. KPS) GO TO 70
  114.             T = Z(KPS)
  115.             Z(KPS) = Z(KP)
  116.             Z(KP) = T
  117.    70    CONTINUE
  118.          IF (Z(K) .NE. 0.0E0) EK = SIGN(EK,Z(K))
  119.          Z(K) = Z(K) + EK
  120.          CALL SAXPY(K-KS,Z(K),A(1,K),1,Z(1),1)
  121.          IF (KS .EQ. 1) GO TO 80
  122.             IF (Z(K-1) .NE. 0.0E0) EK = SIGN(EK,Z(K-1))
  123.             Z(K-1) = Z(K-1) + EK
  124.             CALL SAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1)
  125.    80    CONTINUE
  126.          IF (KS .EQ. 2) GO TO 100
  127.             IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 90
  128.                S = ABS(A(K,K))/ABS(Z(K))
  129.                CALL SSCAL(N,S,Z,1)
  130.                EK = S*EK
  131.    90       CONTINUE
  132.             IF (A(K,K) .NE. 0.0E0) Z(K) = Z(K)/A(K,K)
  133.             IF (A(K,K) .EQ. 0.0E0) Z(K) = 1.0E0
  134.          GO TO 110
  135.   100    CONTINUE
  136.             AK = A(K,K)/A(K-1,K)
  137.             AKM1 = A(K-1,K-1)/A(K-1,K)
  138.             BK = Z(K)/A(K-1,K)
  139.             BKM1 = Z(K-1)/A(K-1,K)
  140.             DENOM = AK*AKM1 - 1.0E0
  141.             Z(K) = (AKM1*BK - BKM1)/DENOM
  142.             Z(K-1) = (AK*BKM1 - BK)/DENOM
  143.   110    CONTINUE
  144.          K = K - KS
  145.       GO TO 60
  146.   120 CONTINUE
  147.       S = 1.0E0/SASUM(N,Z,1)
  148.       CALL SSCAL(N,S,Z,1)
  149. C
  150. C     SOLVE TRANS(U)*Y = W
  151. C
  152.       K = 1
  153.   130 IF (K .GT. N) GO TO 160
  154.          KS = 1
  155.          IF (KPVT(K) .LT. 0) KS = 2
  156.          IF (K .EQ. 1) GO TO 150
  157.             Z(K) = Z(K) + SDOT(K-1,A(1,K),1,Z(1),1)
  158.             IF (KS .EQ. 2)
  159.      *         Z(K+1) = Z(K+1) + SDOT(K-1,A(1,K+1),1,Z(1),1)
  160.             KP = IABS(KPVT(K))
  161.             IF (KP .EQ. K) GO TO 140
  162.                T = Z(K)
  163.                Z(K) = Z(KP)
  164.                Z(KP) = T
  165.   140       CONTINUE
  166.   150    CONTINUE
  167.          K = K + KS
  168.       GO TO 130
  169.   160 CONTINUE
  170.       S = 1.0E0/SASUM(N,Z,1)
  171.       CALL SSCAL(N,S,Z,1)
  172. C
  173.       YNORM = 1.0E0
  174. C
  175. C     SOLVE U*D*V = Y
  176. C
  177.       K = N
  178.   170 IF (K .EQ. 0) GO TO 230
  179.          KS = 1
  180.          IF (KPVT(K) .LT. 0) KS = 2
  181.          IF (K .EQ. KS) GO TO 190
  182.             KP = IABS(KPVT(K))
  183.             KPS = K + 1 - KS
  184.             IF (KP .EQ. KPS) GO TO 180
  185.                T = Z(KPS)
  186.                Z(KPS) = Z(KP)
  187.                Z(KP) = T
  188.   180       CONTINUE
  189.             CALL SAXPY(K-KS,Z(K),A(1,K),1,Z(1),1)
  190.             IF (KS .EQ. 2) CALL SAXPY(K-KS,Z(K-1),A(1,K-1),1,Z(1),1)
  191.   190    CONTINUE
  192.          IF (KS .EQ. 2) GO TO 210
  193.             IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 200
  194.                S = ABS(A(K,K))/ABS(Z(K))
  195.                CALL SSCAL(N,S,Z,1)
  196.                YNORM = S*YNORM
  197.   200       CONTINUE
  198.             IF (A(K,K) .NE. 0.0E0) Z(K) = Z(K)/A(K,K)
  199.             IF (A(K,K) .EQ. 0.0E0) Z(K) = 1.0E0
  200.          GO TO 220
  201.   210    CONTINUE
  202.             AK = A(K,K)/A(K-1,K)
  203.             AKM1 = A(K-1,K-1)/A(K-1,K)
  204.             BK = Z(K)/A(K-1,K)
  205.             BKM1 = Z(K-1)/A(K-1,K)
  206.             DENOM = AK*AKM1 - 1.0E0
  207.             Z(K) = (AKM1*BK - BKM1)/DENOM
  208.             Z(K-1) = (AK*BKM1 - BK)/DENOM
  209.   220    CONTINUE
  210.          K = K - KS
  211.       GO TO 170
  212.   230 CONTINUE
  213.       S = 1.0E0/SASUM(N,Z,1)
  214.       CALL SSCAL(N,S,Z,1)
  215.       YNORM = S*YNORM
  216. C
  217. C     SOLVE TRANS(U)*Z = V
  218. C
  219.       K = 1
  220.   240 IF (K .GT. N) GO TO 270
  221.          KS = 1
  222.          IF (KPVT(K) .LT. 0) KS = 2
  223.          IF (K .EQ. 1) GO TO 260
  224.             Z(K) = Z(K) + SDOT(K-1,A(1,K),1,Z(1),1)
  225.             IF (KS .EQ. 2)
  226.      *         Z(K+1) = Z(K+1) + SDOT(K-1,A(1,K+1),1,Z(1),1)
  227.             KP = IABS(KPVT(K))
  228.             IF (KP .EQ. K) GO TO 250
  229.                T = Z(K)
  230.                Z(K) = Z(KP)
  231.                Z(KP) = T
  232.   250       CONTINUE
  233.   260    CONTINUE
  234.          K = K + KS
  235.       GO TO 240
  236.   270 CONTINUE
  237. C     MAKE ZNORM = 1.0
  238.       S = 1.0E0/SASUM(N,Z,1)
  239.       CALL SSCAL(N,S,Z,1)
  240.       YNORM = S*YNORM
  241. C
  242.       IF (ANORM .NE. 0.0E0) RCOND = YNORM/ANORM
  243.       IF (ANORM .EQ. 0.0E0) RCOND = 0.0E0
  244.       RETURN
  245.       END
  246.